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ABSTRACT 


Models  of  geophysically  important  properties  of  the  Earth,  such  as  seismic  velocity,  Q  and  density,  can  become 
large  and  complex  when  those  properties  vary  in  three  dimensions  within  the  model.  We  have  developed  a  system  to 
represent  the  distribution  of  seismic  properties  in  the  Earth  that  can  accommodate  a  wide  range  of  local  to  global 
scale  3D  Earth  models  with  spatially  variable  resolution.  A  2D  grid  of  nodes  is  tessellated  using  either  triangles  or 
quadrilaterals  and  a  profile  is  defined  at  each  2D  grid  node  that  extends  from  the  center  of  the  Earth  to  the  surface. 
The  surface  of  the  model  corresponds  with  the  topographic/bathymetric  surface  of  the  Earth,  which  is  referenced  to 
the  surface  of  the  GRS80  ellipsoid.  Each  profile  can  be  separated  into  a  number  of  layers  defined  by  interfaces 
across  which  geophysical  properties  may  be  discontinuous.  Within  the  layers  between  interfaces,  the  vertical 
distribution  of  geophysical  properties  may  be  defined  by  a  number  of  continuous  sublayers,  by  arbitrary  order 
polynomials  or  by  various  types  of  splines.  Layer  thicknesses  can  vary  laterally  and  zero  thickness  layers  and  layer 
pinch-outs  are  accommodated.  The  distribution  of  nodes  is  very  flexible,  allowing  model  resolution  to  vary  over  a 
wide  range  in  three  dimensions.  In  this  paper,  we  present  detailed  descriptions  of  the  software  algorithms  used  to 
construct,  store  and  interpolate  these  models. 
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Objective 

Complex  models  representing  the  three  dimensional  distribution  of  seismological  properties  within  extensive  areas 
of  the  Earth  are  becoming  ever  more  common  in  seismological  endeavors.  Notable  examples  include  the  global 
Crust  2.0  model  (Laske  et.  al.,  2008),  the  broad  regional  (Western  Eurasia  and  North  Africa  (WENA)  model 
(Pasyanos  et  al.  2004),  and  WINPAK3D  model  (India-Pakistan;  Reiter  et  al,  2001).  The  prevalence  of  such  models 
has  come  about  in  part  as  a  result  of  the  dramatic  improvement  in  computational  capabilities  in  recent  years. 
Computer  systems  able  to  handle  large  (gigabyte),  geographically  extensive  models  that  represent  properties  from 
deep  in  the  Earth  to  its  surface  are  now  commonly  available. 

Unfortunately,  the  level  of  sophistication  of  model  representation  has  not  kept  pace  with  the  models  themselves. 
Many  models  are  still  based  on  regularly  spaced  grids  that  lead  to  increasing  densities  of  nodes  towards  the  poles 
and  deep  within  the  Earth.  These  simplistic  representations  are  employed  not  because  they  match  the  structure  of  the 
Earth,  but  because  they  are  easy  to  implement  and  use.  However,  these  representations  are  inherently  inefficient  and 
they  are  particularly  ill-suited  to  global  scale  models.  To  develop  global  models  that  can  be  accessed  and  used 
quickly,  a  model  representation  is  needed  that  is  designed  specifically  to  fit  the  real  ellipsoidal,  laterally  and  radially 
heterogeneous  structure  of  the  Earth.  In  this  paper  we  describe  an  Earth  model  representation  that  enjoys  global 
coverage  in  both  radial  and  geographic  dimensions,  and  variable  resolution  that  reflects  the  state  of  knowledge  of 
underlying  parameters  rather  than  ease  of  interpolation.  As  described  in  detail  below,  our  representation  involves  a 
two-dimensional  tessellation  of  triangles  on  the  surface  of  a  unit  sphere,  with  radial  profdes  of  Earth  properties 
defined  at  points  on  the  sphere  where  three  or  more  triangles  intersect.  We  also  describe  the  types  of  information 
that  must  be  stored  in  order  to  completely  represent  the  model  and  methods  for  interpolating  information  from  the 
model. 

RESEARCH  ACCOMPLISHED 
Coordinate  System 

Before  we  can  describe  the  model  representation,  we  first  define  the  Earth- 
centered  coordinate  system  used  to  define  the  positions  of  points  relative  to 
the  Earth.  Points  are  defined  by  a  unit  vector,  x  =  {x0,xh  x2 },  with  its 
origin  at  the  center  of  the  Earth,  and  radius,  r,  measured  in  km  from  the 
center  of  the  Earth.  As  illustrated  in  Figure  1,  this  coordinate  system  is 
oriented  such  that  Xo  points  from  the  center  of  the  Earth  toward  the  point  on 
the  surface  with  latitude  and  longitude  0°,  0°;  Xj  points  toward  latitude, 
longitude  0°,  90°  and  x2  points  toward  the  north  pole. 


Figure  1.  Earth  centered  Cartesian 

Earth  Model  Representation  coordinate  system. 

The  Earth  model  representation  we  describe  here  consists  of  a  two-dimensional,  multi-level  tessellation  of  triangles 
which  completely  covers  the  surface  of  a  unit  sphere  without  gaps  or  overlaps.  The  highest  level  of  the  tessellation 
consists  of  large  triangles  which  are  subdivided  into  ever  smaller  triangles  at  lower  levels  of  the  tessellation.  At  the 
vertices  of  the  triangles,  which  all  reside  at  the  surface  of  the  unit  sphere.  Earth  properties  are  defined  along  radial 
profiles  that  extend  from  the  center  of  the  Earth  to  its  surface.  We  first  describe  how  the  tessellation  is  constructed, 
and  then  how  to  find  an  estimate  of  the  value  of  an  Earth  property  at  an  arbitrary  point  in  the  Earth  by  interpolating 
values  stored  in  the  model. 


Definitions 

To  facilitate  description  and  manipulation  of  the  tessellation,  the  following  terms  are  defined: 
•  If  vertex  V  is  a  corner  of  triangle  T,  then  V  is  a  member  of  T  and  T  contains  V. 
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•  Given  vertex  Vt  that  is  a  member  of  triangle  T,  then  Vj  is  the  next  vertex  in  T  if  one  arrives  at  Vj  by 
traversing  the  edge  of  T  that  leaves  V,-  in  a  clockwise  direction  as  viewed  from  outside  the  unit  sphere. 

•  Triangle  T  has  edges,  Eh  i=  1,3,  where  E,  connects  vertices  Vj  and  Fa-,  where  Vj  is  the  next  vertex  in  T  after 
Vi  and  Ft  is  the  next  vertex  in  T  after  Vj.  Note  that  given  a  vertex  Vh  E;  is  the  edge  of  T  that  does  not  contain 

V, 

•  Triangle  Nt  is  the  fth  neighbor  of  triangle  T  if  N-,  is  the  triangle  on  the  other  side  of  edge  Et. 

Each  vertex  should  maintain  a  triangle  membership  list ,  i.e.,  a  list  of  the  triangles  of  which  it  is  a  member. 

Each  triangle  T  should  maintain  the  following  information: 

1 .  References  to  the  three  vertices  which  are  members  of  T 

2.  References  to  its  three  neighbors 

3.  References  to  three  new  vertices.  These  should  default  to  null,  but  may  be  set  to  instantiated  vertices  if  T  is 
subdivided,  as  will  be  described  shortly. 

4.  A  reference  to  a  descendent  triangle,  if  there  are  any.  If  T  is  subdivided,  as  described  later,  the  reference 
can  refer  to  any  of  the  triangles  that  are  descendents  of  T. 

5.  A  Boolean  flag  indicating  whether  or  not  T  is  ‘marked.’  All  triangles  start  out  ‘unmarked’. 


Tessellation  Construction 

Construction  of  each  level  in  the  tessellation  involves 

1 .  Identifying  the  vertices  involved  in  the  given  level  of  the 
tessellation. 

2.  Connecting  the  vertices  together  to  define  triangles  that 
completely  span  the  surface  of  the  unit  sphere. 

3.  For  each  vertex,  populating  its  triangle  membership  list. 

4.  For  each  triangle,  establishing  which  triangles  are  its  neighbors. 

The  first  level  of  the  tessellation  consists  of  a  regular  icosahedron,  which 
is  a  convex  regular  polyhedron  composed  of  twenty  congruent, 
equilateral  triangles  that  meet  at  twelve  vertices  (Figure  2).  Using  the 
coordinate  system  defined  earlier,  the  coordinates  of  the  twelve  vertices 
of  a  regular  icosahedron  are  given  in  Table  1 . 


Figure  2.  A  regular  icosahedron 
consisting  of  12  vertices 
and  20  congruent 
equilateral  triangles. 


Table  1  -  Vertices  of  a  regular  icosahedron. 


Vertex 

Vo 

Vi 

V2 

0 

0.000000000000000 

0.000000000000000 

1.000000000000000 

1 

0.894427190999916 

-0.000000000000000 

0.447213595499958 

2 

0.276393202250021 

0.850650808352040 

0.447213595499958 

3 

-0.723606797749979 

0.525731112119134 

0.447213595499958 

4 

-0.723606797749979 

-0.525731112119134 

0.447213595499958 

5 

0.276393202250021 

-0.850650808352040 

0.447213595499958 

6 

0.723606797749979 

-0.525731112119134 

-0.447213595499958 

7 

0 . 723606797749979 

0.525731112119134 

-0.447213595499958 

8 

-0.276393202250021 

0.850650808352040 

-0.447213595499958 

9 

-0.894427190999916 

0.000000000000000 

-0.447213595499958 

10 

-0.276393202250021 

-0.850650808352040 

-0.447213595499958 

11 

0.000000000000000 

0.000000000000000 

-1.000000000000000 

These  vertices  are  connected  together  to  form  20  congruent,  equilateral  triangles  as  shown  in  Figure  3.  The  20 
triangles  are  numbered  from  0  to  19  and  the  order  in  which  the  vertices  are  connected  together  is  given  in  Table  2. 
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Table  2  -  Connectivity  of  initial  tessellation. 


Triangle 

Vertex  0 

Vertex  1 

Vertex  2 

0 

0 

2 

1 

1 

0 

3 

2 

2 

0 

4 

3 

3 

0 

5 

4 

4 

0 

1 

5 

5 

2 

7 

1 

6 

3 

8 

2 

7 

4 

9 

3 

8 

5 

10 

4 

9 

1 

6 

5 

10 

1 

7 

6 

11 

2 

8 

7 

12 

3 

9 

8 

13 

4 

10 

9 

14 

5 

6 

10 

15 

7 

11 

6 

16 

8 

11 

7 

17 

9 

11 

8 

18 

10 

11 

9 

19 

6 

11 

10 

Figure  3.  Icosahedral  net  illustrating  one  possible  node 
and  triangle  numbering  scheme. 


As  each  triangle  listed  in  Table  2  is  instantiated,  a  reference  to  the  triangle  is  added  to  each  of  its  vertex’s  triangle 
membership  list. 

Triangle  Neighbor  Identification 

Once  a  set  of  vertices  have  been  combined  into  a  set  of  triangles  at  a  given  level  of  the  tessellation,  and  the  vertices 
have  updated  their  triangle  membership  lists,  it  is  necessary  to  discover,  for  each  triangle  T  at  the  current  level  of  the 
tessellation,  references  to  the  triangles  Nh  i=  1,3,  that  are  the 
neighbors  of  T.  This  is  accomplished  by  performing  the 
following  operations  for  each  vertex  V ,  in  T  (see  Figure  4): 

1. 


2. 

3. 


4. 


Identify  the  other  two  vertices  in  T,  V/  and  Vk,  such  that 
Vj  is  the  next  vertex  after  V,-  and  Vk  is  the  next  vertex 
after  Vk. 

Mark  each  triangle  in  V/s  triangle  membership  list  that 
resides  at  the  same  level  of  the  tessellation  as  T. 

Search  through  the  triangles  in  Vk’s  triangle 
membership  list.  Exactly  two  triangles  will  be  marked: 
the  original  triangle  T,  and  the  neighbor  of  T  that  lies 
across  the  edge  Eh  thereby  allowing  IV,-  to  be  identified. 
Unmark  all  the  triangles  which  were  marked  in  step  2. 


® 

Figure  4  -  Algorithm  for  identifying  triangle 
T  s  fth  neighbor,  Nt.  Triangles  containing  a 
black  dot  are  triangles  that  have  bee  marked. 


Descendent  Level  Construction 

From  a  given  level  of  the  tessellation,  a  lower  level  can  be  constructed  by  subdividing  some  or  all  of  the  triangles  of 
the  given  level  into  smaller  triangles.  Variable  resolution  is  achieved  by  only  subdividing  a  subset  of  the  triangles  at 
the  higher  level  of  the  tessellation.  Construction  of  a  new,  lower  level  from  a  given  higher  level  is  a  three-step 
process. 

Step  1 

Visit  each  triangle  T  at  the  lowest  level  of  the  tessellation.  If  it  is  determined  that  T  should  be  subdivided,  then  visit 
each  edge  Eh  of  T.  If  there  is  not  already  a  new  vertex  positioned  in  the  middle  of  E,.  then 

1 .  Instantiate  a  new  vertex,  V„,  in  the  middle  of  E„ 
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2.  Add  a  reference  to  V„  to  Fs  list  of  new  vertices 

3.  Add  a  reference  to  V„  to  Fs  fth  neighbor,  N,’  s  list  of  new  vertices. 

After  all  the  triangles  at  the  current  level  of  the  tessellation  have  been  visited,  each  triangle  will  have  had  0,  1,  2  or  3 
new  vertices  added  to  its  list  of  new  vertices. 


Step  2 


Each  triangle  T  at  the  current  level  of  the  tessellation  is  visited 
a  second  time  and  the  number  of  new  vertices  added  to  its 
edges  is  identified.  If  three  new  vertices  were  added,  then  T  is 
divided  into  four  new  triangles,  as  illustrated  in  Figure  5a.  If 
two  new  vertices  were  added,  then  there  are  two  distinct  ways 
that  T  can  be  divided  into  three  new  triangles,  as  illustrated  in 
Figure  5c  and  5d.  The  choice  can  be  made  arbitrarily.  If  a 
single  new  vertex  was  added,  then  T  is  divided  into  two  new 
triangles  as  illustrated  in  Figure  5b.  If  no  new  vertices  were 
added  to  T,  then  a  new  triangle  is  instantiated  that  has  the 
same  three  vertices  as  T.  It  is  necessary  to  create  this  new 
triangle,  even  though  it  has  the  same  vertices  as  T,  because  the 
new  triangle  will  ultimately  reside  on  a  different  level  of  the 
tessellation  and  have  different  neighbors  than  T. 


One  of  the  newly  created  triangles  must  be  identified  as  the 
descendent  of  T.  If  T  was  divided  into  three  or  four  new 
triangles,  then  the  center  one  should  become  Fs  descendent. 
Otherwise  there  is  either  no  choice  to  make,  or  the  choice  is 
arbitrary. 


Figure  5.  An  illustration  of  some  of  the  ways  in 
which  a  triangle  can  be  subdivided 
into  smaller  triangles,  a)  three  new 
vertices,  b)  one  new  vertex,  c  and  d) 
two  new  vertices. 


Once  all  the  new  triangles  have  been  created  and  each  old  triangle’s  descendent  identified,  then  all  the  new  triangles 
are  added  to  a  new  level  of  the  tessellation.  The  third  and  final  step  of  the  creation  of  the  new  level  of  the 
tessellation  is  to  determine  each  triangle’s  three  neighbors  as  described  above. 


Property  Profiles  at  Vertices 


At  each  vertex  of  the  tessellation,  a  radial  property  profile  is  defined  which  is  comprised  of  some  number  of  major 
layers.  All  property  values  within  a  major  layer  must  be  radially  continuous  but  property  values  may  be 
discontinuous  across  major  layer  boundaries.  Examples  of  major  property  discontinuities  within  the  Earth  are 
sedimentary  layer  boundaries  in  the  upper  crust;  boundaries  between  the  upper,  middle,  and  lower  crust;  the  Moho, 
the  410  km  discontinuity,  the  660  km  discontinuity,  the  core-mantle  boundary  (CMB),  and  the  inner  core  boundary 
(ICB).  The  radial  property  distribution  within  a  major  layer  can  be  constant,  or  it  can  vary.  Currently,  variations  can 
be  defined  by  a  number  of  sublayers,  by  a  polynomial  in  radius,  or  by  a  cubic  spline.  Support  for  other  distributions 
could  easily  be  added. 

Property  values  within  each  layer,  and  the  radii  of  the  layer  boundaries,  may  be  different  along  profiles  at  adjacent 
tessellation  vertices,  and  layers  may  pinch  out  to  zero  thickness.  Since  linear  interpolation  between  adjacent  profiles 
is  employed,  as  described  next,  property  values  and  interface  radii  will  be  laterally  continuous. 

Model  Interpolation 

Determining  estimates  of  the  property  values  and  layer  radii  at  an  arbitrary  position  Ain  the  Earth  consists  of  the 
following  steps: 

1.  Finding  the  triangle  T  in  the  tessellation  such  that  Aresides  in  T  and  Thas  no  descendents. 

2.  Determining  the  interpolation  coefficients  in  geographic  dimensions. 

3.  Computing  the  major  layer  radii  along  an  interpolated  profile  at  the  geographic  position  of  A  using  the 
interpolation  coefficients  computed  in  step  2. 
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4.  Identifying  the  layer  in  the  interpolated  profde  in  which  the  radius  ofXresides. 

5.  For  each  vertex  Vt  of  T.  interpolating  the  value  of  the  desired  property  at  the  radius  ofX,  but  constrained  to 
be  in  the  layer  identified  in  step  4. 

6.  Computing  the  interpolated  property  values  at  X  from  the  three  property  values  calculated  in  step  5. 

These  steps  are  now  described  in  more  detail. 


Triangle  Walking  Algorithm 

To  find  the  triangle  in  which  Xresides,  a  triangle  walking  algorithm  is  implemented  (Lawson,  1977).  This  method  is 
the  standard  approach  of  point  searching  a  well-behaved,  convex-everywhere,  2D  mesh.  We  initialize  T  to  be  an 
arbitrary  triangle  at  the  highest  level  of  the  tessellation.  For  each  edge  E ,  in  T  we  compute  the  scalar  triple 
product s,  =  {Ek  x  V  )•  X  ,  where  Vj  and  Vk  are  the  first  and  second  vertices  of  Eh  respectively.  X resides  in  T  if  st  >  0 
for  i  =  1,3  .  If  Si  is  negative  for  any  E ,  in  T,  then  T  is  set  equal  to  the  neighbor  that  resides  on  the  other  side  of  edge  is, 
and  the  search  continues.  When  T  is  identified  such  that  st  >  0  for  i  =  1,3 ,  T  is  checked  to  see  if  it  has  a  descendent. 
If  it  does,  then  T  is  set  equal  to  the  descendent  and  the  search  continues.  The  search  ends  when  A  resides  in  T  and  T 
has  no  descendents. 

The  interpolation  coefficients  to  be  applied  to  the  three  vertices  Vt,  1=1,3  are 
1  =  1,3  (1) 

where 

3 

S  =  y\  (2) 

Next,  we  identify  the  3  property  profiles  at  vertices  Vh  and  instantiate  an  interpolated  profile,  /,  at  the  geographic 
position  of  point  X.  Then  we  compute  interpolated  values  for  the  radii  of  all  the  major  layer  boundaries  in  /  using 

3 


where  rj  and  ry  are  the  radii  of  the  f  th  major  layer  boundaries  in  /  and  P„  respectively.  We  then  identify  m,  the  index 
of  the  major  layer  in  I  in  which  the  radius  of  X  resides  using  bisection.  Next,  we  interpolate  property  values  v,  at  the 
radius  of  X  in  profiles  P,  using  whatever  interpolation  algorithms  are  appropriate  for  the  property  distributions  in 
layer  m  of  profiles  P,.  The  interpolations  should  be  constrained  to  layer  m,  i.e.,  if  the  radius  of  A  is  greater  than  the 
radius  of  the  top  of  layer  m  in  P„  or  less  than  the  radius  of  the  bottom  of  layer  m  in  P,.  then  v,  is  set  to  the  value  of 
the  desired  property  at  the  top  or  bottom  or  layer  m  in  Ph  as  appropriate.  Finally,  we  compute  the  interpolated  value 
of  the  property  at  X 

3 
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Example  Model 

To  assess  some  of  the  advantages  and  disadvantages  of  the  unstructured  tessellations  described  here  relative  to  more 
commonly  encountered  structured  grids,  we  have  constructed  an  unstructured  version  of  the  Crust  2.0  Model  of 
Laske,  et.  al.  (2008),  which  is  represented  with  uniform,  2°  latitude  and  longitude  spacing.  We  built  the  unstructured 
tessellation  in  the  manner  described  above,  with  the  criterion  that  triangles  were  only  subdivided  if  they: 

•  contained  at  least  two  of  the  original  nodes  from  the  Crust  2.0  Model; 

•  the  geophysical  parameter  values  at  any  two  of  the  original  nodes  were  different;  and 

•  the  triangle  edges  were  larger  than  2°. 

Our  triangles  are  no  smaller  than  2°  and  at  low  latitudes  they  either  contain  no  original  nodes,  or  all  of  the  original 
nodes  that  they  do  contain  are  identical  to  each  other.  At  high  latitudes,  where  the  Crust  2.0  Model  is  very  highly 
sampled,  many  triangles  do  contain  original  nodes  that  are  not  identical  in  order  to  avoid  violating  the  2°  minimum 
triangle  size  constraint.  Figure  6  shows  an  orthographic  projection  of  the  globe,  centered  on  the  boundary  between 
the  Pacific  Ocean  and  South  America,  illustrating  the  economies  achieved  over  the  oceans  where  the  crustal 
structure  is  constant  over  broad  areas. 


Figure  6  -  Orthographic  projection  of  crustal  thickness  in  the  Eastern  Pacific  and  South  America. 


Figure  7  is  an  azimuthal  equal  angle  projection  of  the  area  near  the  pole,  illustrating  that  the  triangular  cells  maintain 
approximately  constant  area,  even  at  high  latitudes,  while  the  cells  of  the  regularly  spaced,  2°  x  2°  grid  are 
characterized  by  very  high  aspect  ratios. 
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Figure  7.  Comparison  of  the  unstructured  triangular  tessellation  and  the  regularly  spaced  2°  x  2°  original 
grid  at  high  latitude.  The  North  Pole  is  near  the  left  edge  of  the  image. 


Figure  8  illustrates  that  a  2“  x  2°  rectangular 
box  near  60°  latitude  has  half  the  surface  area 
of  a  similar  box  at  the  equator  and  such  a  box 
near  the  pole  has  less  than  2%  of  the  surface 
area  of  the  equatorial  box. 

Conclusion 

We  have  described  a  method  for  representing  t<r- 
complex  3D  Earth  models  based  on  a  multi-  J 

level  triangular  tessellation  of  a  unit  sphere.  At  o 
each  triangle  vertex,  geophysical  properties  are  nT 
defined  along  a  radial  profile  that  extends  from  < 
the  center  of  the  Earth,  through  the  triangle 
vertex  and  out  to  the  topographic/bathymetric 
surface  of  the  Earth.  Construction  of  the 
tessellation  starts  with  an  icosahedron  and 
proceeds  by  triangular  subdivision.  Variable 
resolution  is  accomplished  by  subdividing 
triangles  into  smaller  triangles  only  where  the 
distribution  of  the  geophysical  properties 
represented  in  the  model  warrants  additional 


Figure  8.  Area  at  the  surface  of  the  Earth  of  a  2°  x  2° 
rectangular  box,  as  a  function  of  latitude. 
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resolution.  Unlike  structured,  latitude-longitude  grids,  resolution  need  not  increase  arbitrarily  at  high  latitudes, 
making  this  method  of  model  construction  particularly  suitable  for  global  scale  models. 

Efficient  interpolation  of  geophysical  properties  at  arbitrary  points  in  the  model  is  accomplished  by  implementation 
of  a  walking  triangle  algorithm  which  starts  at  the  highest  level  of  the  multi-level  tessellation  and  proceeds  down 
through  successively  deeper  levels  until  the  desired  position  is  reached. 
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